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; ABSTRACT 

. Accurate phase space coordinates (three components of position and velocity) of individ- 

Oh' ual halo stars are rapidly becoming available with current and future surveys such as SDSS- 

Q ; SEGUE, RAVE, LAMOST, SkyMapper, HERMES and eventually Gaia. This will enable the 

■ computation of the full 3-dimensional orbits of these stars. Spectral analysis of halo orbits 
' can be used to construct "frequency maps" which provide a compact representation of the 
, 6-dimensional phase space distribution function. A frequency map readily reveals the most 
' important major orbit families in the halo, and the orbital abundances in turn, reflect the shape 

^ I and orientation of the dark matter halo relative to the disk. We apply spectral analysis to 

^ , halo orbits in a series of controlled simulations of disk galaxies. Although the shape of the 

• simulated halo varies with radius, frequency maps of local samples of halo orbits confined 
Q>^ I to the inner halo contain most of the information about the global shape of the halo and its 
T— I ■ major orbit families. Quiescent or adiabatic disk formation results in significant trapping of 
CO ' halo orbits in resonant orbit families (i.e. orbits with commensurable frequencies). If a good 

, estimate of the Galactic potential in the inner halo (within ^ 50 kpc) is available, the appear- 

■ ance of strong, stable resonances in frequency maps of halo orbits will allow us to determine 
' the degree of resonant trapping induced by the disk potential. The locations and strengths 

T-H . of these resonant families are determined both by the global shape of the halo and its dis- 

' tribution function. Identification of such resonances in the Milky Way's stellar halo would 

. !^ I therefore provide evidence of an extended period of adiabatic disk growth. If the Galactic 

. potential is not known exactly, a measure of the diffusion rate of large sample of ~ 10'' halo 

' orbits can help distinguish between the true potential and an incorrect potential. The orbital 

\ spectral analysis methods described in this paper provide a strong complementarity to exist- 
ing methods for constraining the potential of the Milky Way halo and its stellar distribution 
function. 

Key words: Galaxy: evolution. Galaxy: formation. Galaxy: halo. Galaxy: kinematics and 
dynamics. Galaxy: structure; cosmology: dark matter; methods: numerical 



1 INTRODUCTION 

Over the last decade a coherent picture of the formation of the stel- 
lar halo of the Milky Way (hereafter MW) has begun to emerge 
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both from observational surveys of the Galaxy and cosmologically 
motivated simulations. In this "concordance cosmological model" 
galaxies are embedded in dark matter halos and form via the merger 
of protogalactic fragments. The fragments consist of dark matter, 
gas, and the first generations of stars which formed in high density 
peaks in the early Universe. An extended period of "late infall" in- 
cluding the ac cretion and tidal shredding of numerous dwarf satel- 
lite g alaxies ( Bullock & Johnston 2005: De Lucia & Helmil l2008l : 
iFont etal. 2008: Jo hnston etal. 2008; Cooper et al. 2010) has con- 
tinued to build up the halo even at the present time. In con- 
tributing to the growth of disk galaxies, gas accretion competes 
with accretions and mergers of satellites which tend to thicken 
and disrupt the d isk. Recent resolved-star all-sky surveys such 
as Sp SS-SEGUE jYannv et alii2009.) , and RAVE cSteinmetz et al.l 
l2006h have revealed that current observations of the stellar halo 
are broadly consistent with the ACDM galaxy formation paradigm. 
For example the detection in SDSS data of num erous substruc- 
tures in the MW halo in the form of tidal strea ms jNewberg et al.l 
l2002l : lYannv et alJl2003l : lBelokurov et al.ll2006l) , the measurement 
of the de gree of dumpine ss in the distribution of stars in the stel- 
lar halo teell et al] |2008^). and the discovery of numero us ultra- 
faint dwarf spheroidal galax ies iWillman et al. 2005; Zuc ker et al.l 
I2OO6I ; iBelokurov et al.ll2007h . have all bolstered evidence that the 
Milky Way's stellar halo was produced in the manner consistent 
with ACDM. 

Although this picture of galaxy formation has been very suc- 
cessful, the availability of phase space coordinates for tens of 
thousands of halo stars from resolved star surveys have revealed 
the inadequacies, and inherent degeneracies, of current analysis 
tools. For instance phase space data are currently compared with 
simulations using simple measures such as orbital eccentricity 



("Sales et al."20 0^, not only yield degenerate results jWilson et al.l 
[2OIO; Dierick x et al. I I2OIOI) , but depend on assumptions regard- 
ing the potential of the Galaxy. Although gas accretion prob- 
ably dominates in disk galaxies it is unclear from the current 
data whether t he MW's inner stellar halo, and thick disk form ed 
purely in situ /Schonrich & Binnevll2009l ; lLoeb man et al."2010!) or 
whether a pre-existing thin disk was heated Ivillalobos & H elmil 
I2OO8I ; iKazantzidis et al. 2008) or if they were entirely created b y 
satellite accretion jAbadi et al.<2003l ; lYoachim & Dalcantoij2008h . 
There is also evidence that the inner stellar halo may be smoother 
(less rms density v ariation) and more metal rich than the outer halo 
I Helmi et ai] |201ll) . In addition, a recent analysis of the full space 
motions of ~ 17, 000 halo star s front the SPSS-SEGU E calibra- 
tion sample dCarollo et al.]|200l 120101 ; iBeers et al.ll2oT ij) suggests 
that the Galactic halo consists of two overlapping components: an 
inner halo, which is rotating in the same direction as the disk and 
an outer halo, with a small retrograde motion. In addition the two 
components have di fferent density profiles, stellar orbits and metal- 
licites (however see !Schoenrich et alj201Ct for a different view.). 

Current analysis tools rely on comparing the average distribu- 
tions of stars in the halo with those arising from simulations are not 



powerful enough to compare the full 6-dimensional phase space 
distribution functions (DF) of 10* — 10'^ halo stars and the infor- 
mation on their metallicities and ages that will bec ome available 
following the launch of Gaia jPerrvman et alj|200ll) . It is impossi- 
ble for current analyses to rule out the possibility that a significant 
fraction of nearby halo stars are not formed in our own potential, 
and have been heated up into a very thick spheroidal inner halo. 
Addressing questions of this kind are of key importance when try- 
ing to constrain galaxy formation models, and in uncovering the 
formation history of the Milky Way, and can only be achieved with 
new techniques that can quantitatively compare the self-consistent 
chemo-dynamical and phase-space distribution function of the stel- 
lar halo with those from models. 

In this paper we show that frequency analysis tools can be ap- 
plied to orbits of halo stars to uncover the phase space distribution 
function of the entire stellar halo. These methods can be applied 
to both the real DF of MW halo orbits and to those from A'^-body 
simulations, permitting detailed quantitative comparisons of the or- 
bit families in the distribution functions. We also show that the halo 
orbit DF (represented by a frequency map) reflects the global shape 
of the dark matter halo and its orientation relative to the disk. 

One of the robust predictions of the ACDM cosmologi- 
cal paradigm is t hat the dark matter halos o f galaxies like the 
MW are tr iaxial (Dubinski & Carlberg 1991; Jing & Sutol l2002l ; 
iBailin & Ste inmetz 2005; AUgood et al. 2006). Orbits of particles 
in triaxial potentials are very different from those in spherical and 
axisymmetric potentials. For instance none of the orbits conserve 
any component of angular momentum. The fraction of such "tri- 
axial orbits" in a potential is a strong indicator of the global shape 
of the potential. If halos are elongated in ways that reflect their ac- 
cretion history and their orientation relative to large scale filaments 
jHelmi et alj|201 ll) . determining the shape and the orbital popula- 
tions could give us important clues to the formation of the MW. 

Recent simulations have shown that when baryons cool 
and condense at the centers of triaxial dark matter halos, 
the ha los become more a xisymm e tric especially a t small 
radii jPu binski & Carlberg 'l99l'; 'Kazantzidis et al. 2004 
iDebattista et al. 2008; Tissera et al. 2010; Kaz antzidis et al 201ij) . 
Two recent studies (Debattista et al. 20081: IValluri et alj I2OI0I) 
analyzed the orbital properties of halo particles in a series of 
TV-body simulations in which baryonic comp onents were grown 
adiabatically in triaxial and prolate halos. IValluri et alj i2010l 
hereafter VIO) showed that although the inner one-third of the 
halo becomes nearly oblate following the growth of a baryonic 
component, only a fraction of the orbits change their true orbital 
characteristics. When the baryonic component is a disk galaxy, 
most of the halo orbits retain a memory of their orbital actions. 
Since the actions are adiabatic invariants, even nearly oblate halos 
can have a significant fraction of box orbits and long-axis tubes. 

This ability of orbits to retain a memory of their initial con- 
ditions has been previously noted in several numerical studies: 
when the potential is changing fairly violently e.g. during a ma- 
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J or merger jValluri et alj|2007ll or duri ng the tidal disruption of a 
satellite in the field of a larger galaxy jHelmi & de Zeeuwll200d) . 
stellar orbits largely conserve their integrals of motion. In the 
latter case Helmi & de Z eeuw ( 200G) showed that the orbital ac- 
tions can be used to recover the relics of the tidal disruption of a 
dwarf satellite by the MW potential. Recently a number of studies 
jMcMillan & BinnevI |2008| ; lOomez & Hdmil I2OI0I : lOomez et al.1 
I2OIOI) have shown that the orbital frequencies (radial oscillation 
or epicyclic frequency) and Jl^ (tangential oscillation or rotation 
frequency) of stars belonging to a single satellite are tightly cor- 
related and that this correlation remains strong long after the rem- 
nant has become completely well mixed in configuration (physical) 
space. Consequently, correlations between orbital frequencies can 
be used to identify stars that belong to individual accretion events, 
and simultaneously determine the true galac tic potential and the 
time since each satellite galaxy was disrupted jMcMillan & BinnevI 
[2OO8 : Gomez & Helmi 2010). 

Following t he method outlin ed by ICarpintero & AguilaJ 

( Il998l) , VIO (also lDeibele"tai]|20 111) used the oscillation frequen- 
cies of orbits of particles in simulated halos to classify orbits into 
major families. VIO also used orbital frequencies to quantify the 
shapes of orbits (and to relate the orbital shapes to the shapes of 
the halos), to identify orbits which are chaotic, and to identify im- 
portant resonances (regions of phase space occupied by orbits with 
commensurable orbital frequencies). In this paper we show that the 
relative distribution of orbital types and the identification of impor- 
tant resonances populated by halo orbits strongly reflects the orien- 
tation of the triaxial halo relative to the galactic disk. We will show 
that the discovery of large numbers of halo stellar orbits trapped in 
resonances could put constraints on the form of the Galactic poten- 
tial and the DF of the stellar halo. 

Why are resonances important? When a time-dependent force 
acts on a system like a galaxy, long-lived resonant interactions play 
an important role in the evolution of the system and leave imprints 
in its phase space structure. The identification of stars trapped in 
resonances can put constraints on the potential and its components 
(e.g. the bar and spiral arms). Resonant interactions between in- 
dividual orbits and a changing potential influence time-dependent 
(secular) evolution: e.g. it can cause stars to "levitate" to form a 
thick disk (Sridhar & Touma 1996), it may result in the form ation 
of polar rings or counter-rotating disks dTremaine & Yil2000l) . and 
it can cause resonant shocking (or torq uing) of stars in s atellites as 
they are disrupted in dark matter halos ( IChoi et al.l2009l) . "Capture 
into resonance" is studied extensively in the context of planetary 
systems, where a migrating planet can capture planets or planetesi- 
mals into mean motion resonances jMalhotrj 19931 : 1 Yu & Tremaind 
I2OOII ). In the planetary dynamics literature it has been shown that 
resonant trapping of planetesimals occurs in slowly varying poten- 
tials, but can be prevented when the drift or migration rate is suffi- 
ciently high (i.e. non adiabatic) (e.g. ,Ouillen 2006) . The identifica- 
tion of a significant fraction of resonantly trapped halo orbits could 



therefore provide clues to the way in which the halo potential has 
changed over time. 

Laskar's frequency analysis method is particularly good at 
identifying resonances ( Ro butel & L askar 2001). We show that 
when the method is applied to orbits in a self-consistent distribu- 
tion function, the method permits easy identification of the major 
orbit families and an assessment of the relative importance of each 
family to the phase space DF. 

This paper is organized as follows: Section |2] describes the 
simulations analyzed in this paper and briefly describes the fre- 
quency analysis of orbits. Section|3]presents the results of our anal- 
ysis of adiabatic simulations of isolated galactic potentials consist- 
ing of both axisymmetric and triaxial dark matter halos with parti- 
cle di sks of various orientations (taken mostly from jPebattista et al] 
I2OO8I hereafter DOS). We also discuss the effects of disk galaxies 
which form self-consistently from hot halo gas in a spherical or pro- 
late halo. These simulations include the hydrodynam ical effects of 
gas cooling, star formatio n and supernova feedback dStinson et al.l 
l2006l : lRoskar et alj|2008l) . In Section|4]we show how the mean or- 
bital diffusion rates of a large ensemble of orbits selected from a 
distribution function can be used to assess how much it deviates 
from self-consistent equilibrium and discuss how this measure is 
affected when the assumed Galactic potential is incorrect. In Sec- 
tion|5]we summarize our results, and discuss their implications for 
future large data sets which will obtain the 6 dimensional orbits of 
stars in the MW halo. 



2 SIMULATIONS AND NUMERICAL METHODS 
2.1 Simulations 

We analyzed two types of controlled simulations (a) A'^-body sim- 
ulations in which an exponential, thin stellar disk (consisting of 
collisionless particles) was grown adiabatically inside an isolated 
halo; (b) A'^-body-l-SPH hydrodynamical simulations of the forma- 
tion of stellar disks from initially hot gas distributed inside a spher- 
ical or prolate halo and allowed to cool and form a disk of gas in 
which stars form. In general we assume that the potential of the 
galaxy model is perfectly known, and we characterize the distri- 
bution function, its orbital properties, and its dependence on the 
radial variation of the shape of the halo and its orientation rela- 
tive to that of the disk. In practice, the potential and the DF of the 
halo need to be determined simultaneously, or (in the absence of 
adequate kinematical constraints) the potential alone will be de- 
termined, within some radius. We therefore ran a few additional 
simulations, designed to determine if it is possible to constrain the 
potential from the orbital properties of an ensemble of halo orbits. 

In the controlled simulations presented in this paper, initially 
spherical isotropic NFW (Navarro et al. 1996) halos we re gener- 
ated via Eddington's formula ( Binnev & Tremaind l2008l . § 4.3.1) 
with each halo composed of two mass species arranged on shells. 
The inner shell has less massive particles than the outer one, which 
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Run 


»'200 


A/200 


Mb 


h 




Run Description 


Reference 


Name 


[kpc] 


[IQI^Mq] 


[1O"M0] 




[Gyr] 






SNFWD 


85 


0.66 


0.66 


0.1 


5 


Spherical halo+stellar disk 


This paper 


SAl 


215 


4.5 


1.75 


0.039 


5 


Triaxial halo+ short-axis stellar disk 


DOS 


LAI 


215 


4.5 


1.75 


0.039 


5 


Triaxial halo+long-axis stellar disk 


DOS 


lAl 


215 


4.5 


1.75 


0.039 


5 


Triaxial halo+intermediate (y)-axis stellar disk 


DOS 


TAl 


215 


4.5 


1.75 


0.039 


5 


Triaxial halo+tilted stellar disk 


DOS 


SA2 


215 


4.5 


0.52 


0.012 


1.5 


Triaxial halo+barred stellar disk 


DOS (there labelled BAl) 



Table 1. The collisionless simulations in this paper. Mf, is the mass in baryonic disk and is the baryonic mass fraction, tg is the time during which the 
baryonic disk is grown. In all the models the exponential disk has a radial scale length of 3 kpc. The last column contains the references where the simulation 
was first reported. 



allows for higher mass resolution at small radii. Most of the dark 
matter particles in the inner part of the halo have masses of IO'^Mq . 
Prolate and triaxial halos (consisting of 4 x 10*' particles) were gen- 
erated via mergers of the spherical NFW halos (see DOS for details). 

In the simulations, disks of particles were grown (starting from 
nearly zero initial mass) adiabatically and linearly on a timescale 
tg inside a dark matter halo (see Tab. [T] for details of parameters 
of the simulations). Disks were grown in a spherical halo (model 
SNFWD) and in triaxial halos with the disk plane oriented in var- 
ious ways relative to the principle axes of the halo: (a) perpendic- 
ular to the short axis (model SAl), (b) perpendicular to the long 
axis (model LAI), (c) perpendicular to the intermediate axis (model 
lAl), and (d) tilted at an angle of 30° to the x — y plane of the triax- 
ial halo, by rotating it about the y axis (model TAl IB DOS showed 
that in all of these models, the shape of the halo within the inner 
1/3 of the virial radius becomes nearly (but not exactly) oblate fol- 
lowing the growth of the disk, with the short-axis of the oblate part 
of the halo co-aligned with the spin axis of the disk. In most of the 
simulations the disk particles remain stationary, hence the disk is 
rigid throughout. In one case (model SA2), a "light" disk of parti- 
cles was made "live" after the disk had grown to its final mass. This 
disk subsequently formed a bar which per sisted for ^ 10 Gyr be- 
fore d issolving because of the triaxial halo jBerentzen & Shlosmanl 
I2OO6I) . Additional details may be found in DOS. All the collision- 
less simulations were evolved with PKDGRAV, an efficient, multi- 
stepping, parallel tree code (Stadel 2001). Dark matter particles had 
a softening parameter e = 0. 1 kpc, and that of stars was in the range 
e = 60 - 100 pc. 

Two of our controlled simulations contain a baryonic 
(gas-l-star) disk, which forms self-consistently from hot gas in a 
spherical (model SNFWgs) or prolate (model SBgs) halo. The 
baryonic component is 10% of the total mass and initially has 
the same density distribution as the dark matter particles. The 
halo and gas particles are given an initial specific angular mo- 
mentum j, determine by overall cosmological spin parameter A — 
{j/G){\E\/M''^Y''^ = 0.039, which is motivated by cosmolog- 



^ Note that x, j/, z are defined to be the long, intermediate and short axes 
respectively, of the triaxial halo. 



ical A'^-body experiments teullocketal]|200lh . Both the spheri- 
cal halo and the progenitor halos that were merged to produce a 
prolate halo had the same angular momentum parameter A. Each 
component is modeled with 10^ particles, with the dark matter par- 
ticles of mass 10** M0 , and gas particles having an initial mass 
of 10^ M0 . The gas particles are allowed to cool and form stars 
of ty pical ma s s arou nd 3 x lO** Mq following the prescription in 
iStinson et ai] ( I2OO6I) . The net angular momentum allows the gas 
to form a disk as it cools, resulting in a stellar disk as star for- 
mati on occurs. The simulation closely follows that described in 
Rosk ar et all ( |2008|) and is evolved wi th the parallel 7V-body+SPH 
code GASOLINE jWadslev et aD '2004) for 10 Gyr. 

Although controlled simulations are useful for testing the ef- 
fects of disks with different orientations on the halo DF, these mod- 
els are not fully realistic depictions of how disk galaxies form. In 
the current hierarchical structure formation paradigm, disk galaxies 
probably experienced several gas rich mergers, at least in their early 
history, and continue to accrete small satellites today which add to 
the stellar halo. We defer the study of stellar ha los from cosmologi- 
cal simulations drawn from the MUGS project jStinson et alj2O10h 
to a future paper (, Valluri et al.ii201 1,) . 

2.2 Selecting the halo orbit samples 

In each of the simulations we selected 1 — 2 x 10* dark matter par- 
ticles. The particles were either randomly distributed within some 
spherical volume of radius Vg centered on the model's galactic cen- 
ter, or within a region of radius Rs from the location of the "sun" 
(which was assumed to be at S kpc from the Galactic center). Since 
most of the potentials we studied are non-axisymmetric, the az- 
imuthal location of the "sun" in the equatorial plane relative to the 
major axis of the triaxial halo is an additional parameter. Rather 
than choosing a specific (arbitrary) angle in azimuth for the so- 
lar location, we select particles within a "torus" of width Rs and 
with radius S kpc. When we selected subsamples of the torus re- 
gion we found that the results did not depend much on the precise 
azimuthal location of the "sun", so long as the radial region sam- 
pled was 10 kpc in size, or larger. 

Following VIO we studied orbits of halo particles after the 
disk had grown to its final mass and the halo had relaxed to a new 



© 201 1 RAS, MNRAS Q0Q.mi23] 



Orbital spectral analysis of the stellar halo 5 



equilibrium. The orbits were integrated for 50 Gyr in a frozen po- 
tential corresponding to the new equilibrium potential given by the 
full mass distribution of the simulation (dark matter and baryons) 
using an integration scheme based on the PKDGRAV tree. The 
50 Gyr integrations are used only because they yield highly accu- 
rate frequencies for the halo orbits with the longest periods, and 
should not be construed to imply that this is a physically meaning- 
ful time period. (In a similar way, it is possible to integrate orbits of 
MW halo stars in a fixed potential for equally long times in order 
to characterize the nature of their current orbit.) None of the con- 
trolled simulations include a stellar halo, so we assume that the dark 
matter halo orbits can be used to represent the orbits of particles in 
the stellar halo. While this is clearly not an ideal assumption, many 
studies of the stellar halo use prescriptions to "tag" the most tightly 
bou nd dark matter particles in dwarf satellites as "star particles" 
(e.g. lBuUock & Johnstonll 2005l : ICooDer et alJbOld) . Elsewhere we 
analyze a fully cosmological hy drodynamical simu lation of a disk 
galaxy from the MUGS project l lStinson et al.ll2O10h and show that 
in the inner halo, (a region on which we focus) star particles and 
dark matter particles have simi lar orbital properties, justifying this 
assumption jValluri et al J201 ll) . 



2.3 Laskar Frequency Mapping 

In three dimensional Hamiltonian potentials, the phase space struc- 
ture of regular orbits can be described by three actions Ja (a = 
1,2,3) and three angle variables 6a, which constitute a canoni- 
cally conjugate coordinate system. The actions are integrals of mo- 
tion and are conserved along the orbit, while the angles increase 
linearly with time. The angle variables at any time t are given by 
da{t) = 9a{0) + Qat, whcrc the Qa are called "fundamental fre- 
quencies". The actions are adiabatic invariants and consequently 
remain constant as the potential of the system changes adiabati- 
cally. Regular orbits in 3-dimensional potentials can therefore be 
thought of as occupying the surfaces of 3 dimensional tori, with the 
size of each torus characterized by the actions Ja , while the angle 
variables 6a {t) represent the traversal of the orbit over the surface 
of the torus in each dimension. 

Since the angle-action coordinates are related to the classical 
spatial coordinates and momenta via a coordinate transformation, it 
can be shown that traditional space and velocity coordinates can be 
represented by a time series of the form: x{t) = ^71^6''^'=*, and 
similarly for other phase space coordinates e.g. y{t),Vx{t) (where 
the sum is over all terms in the spectrum). A Fourier transform of 
such a time series will yield the spectrum of orbital frequencies ujk 
and associated amplitudes ^fc, that govern the motion of th e or- 
bit teinnev & Spergellll982l.ll984l : lBmnev & Tremaindlioog) . For 
most regular orbits only three of the frequencies in the spectrum (of 
a given orbit) are linearly independent, (i.e. all other frequencies uj^ 
can be expressed as Uk = Ik^i +mkQ2+nk^3, where Ik, mt, rik 
are integers), fii, ^2, ^^3 are nothing other than the "fundamental 



frequencies" described above and the associated 3 orbital actions 
Ji,J2, Ja can be computed from the amplitudes Ak . 

If all orbits in the potential are regular and the DF can be 
written as a continuous function of global actions Ji, J2, Js e.g. 
fiJi, J2, Ja) then the corresponding fundamental frequencies and 
the associated angle variables 61,62, 63 vary in a continuous man- 
ner across phase space. However it is only possible to write the 
DF as a function of global actions for special cases e.g. strictly 
axisymmetric potentials or separable triaxial (Stackel) potentials. 
If analytic global actions are known then one can simply use sur- 
faces of section to map the phase space structure and to identify 
resonant orbits and chaotic orbits (BT08). In the study of three di- 
mensional orbits in realistic potentials analytic integrals of motion 
other than the Hamiltonian are rarely available. If we compute the 
orbital frequency spectrum of an orbit in an arbitrary coordinate 
system e.g. in Cartesian or cylindrical coordinates the resulting or- 
bital frequencies fia; , f2y , fi^ are in no way "f undamental" t o the 
nature of the potential. However as shown by ( lLaskailll990l) any 
canonically conjugate pair of variables can be combined to obtain 
a frequency spectrum, when there are no global actions. These fre- 
quencies are still referred to as "fundamental frequencies" since for 
any regular orbit all the components of the frequency spectrum are 
linear integer combinations of the 3 fundamental e.g. (Qx , fiy , Jlz) 
or (Q,R, Qcj,, fiz), depending on the coordinate system selected. 

iLaskaj ( Il99(]|. [1993) developed a very accurate numeri- 
cal technique "Numerical Analysis of Fundamental Frequen- 
cies" (NAFF) to recover frequencies in completely general po- 
tentials. We use an im plementation of this algorithm due to 
IValluri & MerrittI (1 1998) which was adapted for application to or- 
bits in A'^-body potentials by VIO. 

The NAFF algorithm recovers orbital frequencies from 3 com- 
plex time series consisting of pairs of phase space variables. For tri- 
axial potentials we use a Cartesian coordinate system centered on 
the center of the galaxy and oriented such that x,y,z correspond to 
the major (long), intermediate and minor (short) axes of the poten- 
tial, respectively. In the Cartesian coordinates the Fourier analysis 
is performed on 3 time series of the form fa{t) = ce(t) + ivait) 
(where a — x,y, z). 

For potentials which are axisymmetric or nearly so, most of 
the orbits are tubes which circulate about the symmetry axis. In 
such potentials it is preferable to work in cylindrical polar coor- 
dinates. We transform from the planar coordinates x,y,Vx,Vy to 
plane polar coordinates R, vr, (f), Q, where R — \J x"^ -\- j/^, az- 
imuthal angle = arctan(a:/j/), vr — {xvx + yvy)/R and 
G = xvy — yVx- R,vr are the canonically conjugate radial co- 
ordinate and momentum and hence can be used to define a com- 
plex time series fR{t) = R{t) + ivR{t). However, since </> and Q 
are the angular coordinate and momentum (and not linear coordi- 
nate and momentum) this pair cannot be used to construct the com- 
plex time series used by the freq uency analysis method. Follow- 
ing Papaphilippou & Laskarl ( 1 1996) we use Poincare's symplectic 
polar variables y/2Q cos 4> and \/2Q sin 0, to define the function 
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= V20(cos (f> + i sin (f)) . For motion perpendicular to the equa- 
torial plane we use the complex series fz{t) = z (t) + ivz (t) . 

VIO demonstrate that fundamental frequencies of orbits in 
a self-consistent DF can be used to construct a "frequency map" 
which gives a picture of the phase space structure based on its or- 
bital content. A frequency map of phase space is obtained by plot- 
ting the ratios of fundamental frequencies (e.g. in Cartesian coor- 
dinates: flx/i^z vs. ily/ilz) for a very large number of orbits. 

As we noted above the coordinate system selected to integrate 
the orbits and compute the frequencies is determined more by con- 
venience and convention than any fundamental property of the or- 
bital frequencies obtained. Nonetheless we will see that for triaxial 
systems, the choice of a Cartesian coordinate system aligned such 
that the global principal axes of the model coincide with the coor- 
dinate system, result in different orbit families appearing in distinct 
groups or lines on a frequency map. The use of an appropriate co- 
ordinate system also allows one to identify truly resonant orbits 
families. 

Resonant orbits are regular orbits that have fewer than 3 lin- 
early independent fundamental frequencies which are related via 
integer linear equation such as: lQ,i + mQ.2 + nO.^ = 0, where 
{I, m, n) are small integers. A frequency map can be used to eas- 
ily identify the most important resonant orbit families, since such 
orbits populate straight lines on such a map. The strength (or im- 
portance) of a resonance can be assessed from the number of orbits 
associated with a particular resonance. 

iMerritt & Valluril ([1999*) showed that perfectly resonant or- 
bits in 3-dimensional potentials have two non-zero fundamental 
frequencies and occupy thin two-dimensional surfaces (generally 
multiply connected), in configuration (physical) space. They are 
surrounded by a resonance region consisting of orbits which share 
the oscillation frequencies of the perfectly thin resonant parent, but 
have a third non-zero frequency, which is small but increases as 
the orbit deviates from its resonant parent. Most of these slightly 
non-resonant orbits also appear along the resonance lines in the 
frequency map. Unstable resonances appear as blank lines or blank 
spaces on the frequency map. 

iLaskar et al. also showed that since the frequencies of 

regular orbits can be recovered with very high accuracy, chaotic 
orbits can be easily identified, since their frequencies do not re- 
main constant but drift when computed over two adjacent time 
intervals. They showed that the rate at which the orbits diffuse 
in frequency space is correlated with their degree of stochastic- 
ity. VIO showed that this way of measuring stochasticity was par- 
ticularly useful in Ai'-body potentials (and superior to the better 
known "Lyapunov exponent") since it is able to distinguish be- 
tween diffusion due to micro-ch aos that arises due to discreteness 
noise l lKandrup & Siderisll2003l) and genuinely irregular behavior. 
We refer the reader to VIO for more details. For each orbit we di- 
vide the integration time of 50 Gyr into two consecutive segments 
(ti and t2) and use NAFF to compute the fundamental frequencies 
na{ti),^la{t2)- The "diffusion" rate for each frequency compo- 



nent is then computed as: 

log(A/.)^log| ""^^^)^-^^^^"^^^^ (1) 

We define the diffusion rate for an orbit, log(A/) (logarithm to 
base 10) to be the value associated with the frequency compo- 
nent Qa with the highest amplitude Aa measured over the entire 
time interval (ti + t2)El The larger the value of the diffusion rate 
log(A/), the more chaotic the orbit. We use the diffusion rate to 
distinguish between regular and chaotic orbits, and to distinguish 
weakly chaotic orbits from strongly chaotic orbits. It is important 
to note that for most systems there is a continuous and nearly Gaus- 
sian distribution of diffusion rates (VIO). 

Laskar's method recovers orbital frequencies of regular orbits 
with very high accuracy in ~ 20 — 30 orbital periods making this 
method particularly valuable for studying Galactic halo stars. We 
integrated all orbits for 50 Gyr, but only present results for those 
with orbital periods shorter than 2.5 Gyr (i.e. those that execute 
more than 20 orbital periods). All frequencies in this paper are re- 
ported in units of Gyr~^. 



3 RESULTS OF CONTROLLED SIMULATIONS 
3.1 Nearly oblate axisymmetric halos 

In oblate axisymmetric potentials most orbits are short axis tubes, 
hence orbits are best studied in cylindrical coordinates. We begin 
with the study of the simplest halo DF: a spherical NFW halo in 
which a stellar disk was grown adiabatically (model SNFWD in 
Tab. [TJ. Figure [T] shows histograms of frequencies for 2x10^ 
halo particles with initial rg < 200 kpc, in Cartesian coordinates 
flj;, Sly, Q,z (top panels) and in cylindrical coordinates about the z- 
axis Qh, Qcj, ilz (bottom panels). The initial DF was spherical and 
isotropic and hence the frequency distributions in the three Carte- 
sian directions are identical (dot-dash curves). The growth of an 
axisymmetric disk (rotating about the z-axis) increases the depth 
of the potential and makes the halo flatter in the inner regions. This 
results in an increase in all the orbital frequencies, but axisymme- 
try implies identical distributions for ^lx and Sly. Since the disk 
potential is significantly flatter in the z direction, there is a greater 
increase in Sl^ for particles that lie closer to the center of poten- 
tial (highest values of Qz), accounting for the slight increase in the 
weight of the high frequency tail of the 0,^ distribution. In cylin- 
drical coordinates the frequency Sl,^ (bottom row, middle column) 
describes the motion in the azimuthal direction and is either pos- 
itive or negative depending on whether the orbit rotates counter- 
clockwise or clockwise, about the z-axis. Since the halo was set up 

^ Note: this definition of "diffusion rate" differs slightly from VIO who 
used the value associated with the largest fundamental frequency. The new 
definition was found to more accurately identify chaotic orbits and yields a 
lower rate of misclassification of regulai' orbits as chaotic. 
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Figure 1. Histograms of orbital frequencies of 10^ halo particles with rg < 200 kpc. The vertical axes have arbitrary scales. Top row (L-R): histograms of 
^Ix , , f2z in a spherical isotropic NFW halo (dot-dash curves), and after a thin collisionless stellar disk was grown adiabatically (solid curves); Bottom row 
(L-R): Qn, Q^, and for the case when a disk is present. The distribution in is bimodal because the original halo had a spherical isotropic DF. 



to be nonrotadng, fl^ values are symmetrically distributed about 
zertifl. 

Figure |2] shows frequency maps for the same 2 x 10'* halo 
particles in Cartesian coordinates (left) and in cylindrical coordi- 
nates (right). Each particle is represented by a single point whose 
location is determined by the ratio of the fundamental frequencies 
in Cartesian coordinates (i^x/^z vs. Qy/Q,z) or cylindrical coordi- 
nates (fiz /fifl vs. I Qr). In frequency maps the color of a point 
represents the binding energy of its orbit with blue representing the 
l/3rd most tightly bound particles in the map, red representing the 
l/3rd least bound particles and green representing the intermediate 
energy range. 

Since the growth of the disk makes the originally spherical dis- 
tribution of particles oblate axisymmetric, the DF is entirely popu- 
lated by short-axis tube orbits. In a Cartesian frequency map (Fig- 
ure |2] left panel) such orbits primarily lie along a diagonal line that 
satisfies the condition Q,x/^z ~ Qy/^lz- Most short axis tubes 
are not "resonant" orbits, but in a Cartesian frequency map they all 
appear clustered along a line, because each short-axis tube can be 



^ Note that the frequency fi^ is measured relative to the center of the cyUn- 
drical coordinate system and not relative to a circular orbit as in the case of 
the epicyclic frequency for disk stars. 



viewed as arising from a radial perturbation of a parent "thin shell" 
orbit which is a resonant orbit ( de Zeeuw & Hunter 1990.). 

Truly resonant short-axis tube orbits are those that appear 
clustered along lines in a frequency map in cylindrical coordi- 
nates (Figure [2] right panel). This map shows bisymmetry be- 
tween the right and left halves, reflecting the bimodal distribu- 
tion in (Fig. [T). The frequency map in cylindrical coordi- 
nates shows a striking number of resonances which appear (pri- 
marily) as horizontal lines delineated by the enhanced cluster- 
ing of particles at resonances between the vertical oscillation fre- 
quency flz and radial oscillation frequency Qr. Resonances are 
seen at Qz/ilR = 0.5,0.66,0.75,0.83,1,1.5 (i.e. Uz/nR = 
l/2;2/3,3/4, 5/6, 1/1,3/2), as indicated by labels. Thus the 
growth of a disk in a spherical potential results in halo orbits be- 
coming resonantly trapped at numerous resonances, primarily be- 
tween the radial and vertical oscillation frequencies. 

3.2 Triaxial halos with rigid disks 

We now consider controlled simulations of galactic potentials with 
triaxial halos. The simulations are taken from DOS. We consider 
four different orientations of the disk relative to the same triaxial 
halo (Halo A from DOS). In all four cases the disk potential con- 
sisted of particles, but the particles were held rigid while the halo 
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Figure 2. Frequency maps of 10^ halo orbits in tlie model with a stellar disk grown in a spherical NFW halo. Left: in Cartesian coordinates the map shows 
that most orbits lie along the diagonal "resonance line". This is not a true resonance but represents all orbits associated with the short-axis tube family. Right: 
Frequency map of the same particles in cylindrical coordinates. The bisymmetry about Q^/Q,r = results because the halo has no net rotation. The map 
shows several resonances which appear as horizontal lines: e.g. = 0.5, 0.66, 0.75, 0.83, 1, 1.5. In both panels (and hereafter) particles are color 

coded by binding energy in three energy bins, each containing l/3rd of the particles. 




R [kpc] 

Figure 3. Halo axis ratios b/a (solid curves) and c/a (dashed curves) as a 
function of radius. The green curves show the shape of the initial triaxial 
halo. The black curves are for the triaxial halo -I- short axis disk (model 
SAl), the blue curves are for the triaxial halo -I- long axis disk (model LAI), 
the red curves are for triaxial halo -I- intermediate axis disk (model lAl) and 
the grey curves are for the triaxial halo -I- tilted disk (model TAl). 



was allowed to relax. In all four cases, the inner region of the halo 
became less triaxial (more oblate) with short axis aligned with the 
syrmnetry axis of the disk. We define a, h, c to be the density semi- 



axes of the major, intermediate and short axis respectively of the 
global halo. 

The shapes of the triaxial halos were measured as described 
in DOS. Briefly, we measure the eigenvalues of the unw eighted mo- 
ment of inertia tensor / obtained in bins of particles l lKatjl99ll 
see also eq. 2 from DOS). The ratios of the eigenvalues of the di- 
agonalized moment of inertia tensor (In > T22 > I33) are used 
to calculate the axis ratios b/a and c/a. The shapes were measured 
using an iterative procedure in annular shells of fixed semi-major 
axis width and are differential rather than integrated over all parti- 
cles within a given ellipsoidal radius. 

Figure [3] shows the change in the axes ratios b/a and c/a as 
a function of radius. The green curves show that the initial triaxial 
halo A was very strongly prolate-triaxial (a result of low angular 
momentum mergers). The black curves are for model SAl (triaxial 
halo + short axis disk), the blue curves are for model LAI (triaxial 
halo + long axis disk), the red curves are for model lAl (triaxial 
halo -I- intermediate axis disk), and the grey curves are for model 
TAl (triaxial halo -1- tilted disk). Regardless of the orientation of 
the disk, we see that its growth results in a very significant increase 
in axis ratios b/a (solid curves) and c/a (dashed curves) within the 
inner 30 kpc, for models SAl, LAl, TAl, and a moderate increase 
in oblateness within 50 kpc for model lAlQ. (DOS shows the change 



* Recall that as 6/a — > 1, a model becomes more oblate. 
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in halo shape out to 200 kpc.) Model LAI (blue curves) shows 
the most significant change in shape over the radial range ~ 20 — 
50 kpc. In the inner 15 kpc, the model is more triaxial than it is at 
larger radii, due to the formation of an iimer elongation along the 
y-axis (see Fig.|4}. 

Although all the models become more oblate at the center, 
they still remain triaxial both at small and large radii. This is seen 
in Figure|4]which shows contour plots of the dark matter halo pro- 
jected density (black curves) in two projections for each of the 
models studied. The top row shows density contours of the dark 
matter particles when viewed with disk edge-on as represented by 
the red contours. The bottom row shows the density contours of 
the halo in the plane of the disk, except in model TAl (where the 
disk is inclined to the principle planes). The top panels show that 
in all cases the contours become flattened with the symmetry axis 
co-aligned with the disk's symmetry axis at radii < 15 kpc, but 
retain their original elongation along the x-axis at larger radii. The 
lower panels show that although the halo becomes flattened in the 
inner region, it is not axisymmetric, even in the plane of the disk, 
showing that triaxiality varies with radius. 

In the rest of this section we will show frequency maps of or- 
bits in each of these models. In all cases we show frequency maps 
in Cartesian coordinates for ~ 10* halo orbits selected in two ways. 
First, orbits were randomly selected to have an initial distance from 
the galactic center Vg < 200 kpc. Since these frequency maps rep- 
resent a randomly selected subsample of orbits drawn from the en- 
tire dark matter halo we will refer to these maps as showing the 
"full halo DF". 

Ideally we would like to be able to distinguish between dif- 
ferent halo orientations relative to the disk from the DF of halo 
stars since the orbits of dark matter particles are not directly ob- 
servable. However, halo stars are known to be significantly more 
centrally concentrated than the dark matter (Battaglia et al. 2005). 
Therefore we also present frequency maps of orbits selected in a 
second way, where the orbits satisfy two conditions: (i) their ini- 
tial distance from the "sun" Rs < 10 kpc, (ii) each orbit has an 
apocenter radius rapo < 50 kpc from the galactic center. Restric- 
tion ( i) is motivated by the expectation that Gaia jPerrvman et al.l 
l200lh will obtain 6 phase space coordinates (and in particular the 
most accurate distances and proper motions) for stars within 10 kpc 
of the sun. Restriction (ii) is imposed because we do not anticipate 
that any method of halo shape determination will provide an ac- 
curate measurement of the shape of the halo at distances greater 
than 50 kpc from the Galactic center. By restricting ourselves to 
only those stars with rapo < 50 kpc, we are ensuring that the maps 
only contain orbits which explore the region of the halo where the 
shape is well determined. We will refer to this second set of maps 
as showing the "inner halo DF". Figures [3] & |4]show that the in- 
nermost region of the halo becomes significantly oblate due to the 
presence of the disk, regardless of orientation of the large scale 
halo. 

In the sections that follow we will compare the DFs of the 



four model by showing that the frequency map representations of 
their DFs differ significantly from each other regardless of whether 
the full halo or inner halo is represented. This is remarkable since 
the initial DFs of each halo (prior to the growth of the disk) were 
identical and Figure|4]shows that while the disk alters the inner re- 
gions of the halo it remains triaxial and elongated along its original 
long-axis at large radii. 

It is important to note that in all cases we begin with exactly 
the same set of 10* orbits from the original triaxial halo (whose 
frequency map is plotted in Fig.|5]left), and follow that same set of 
orbits in the different models in which disks were grown adiabati- 
cally. 



3.2.1 SAl: Triaxial halo with short axis disk 

Figure [5] (left) shows the frequency map of a triaxial halo that was 
formed from multiple mergers of spherical NFW halos (This model 
is the original triaxial halo in which all the disks are grown and is 
referred to as "Halo A" in DOS and VIO.) The particles are color 
coded by energy as described above. 

VIO used the relationships between the fundamental frequen- 
cies of orbits in this triaxial potential to characterize them by orbital 
type. They showed that 86% of the particles were on box orbits, 
11% were on long-axis tube orbits, 2% were short-axis tubes and 
1% were chaotic. The fact that box orbits dominate over long-axis 
tubes and short axis tubes in this model strongly reflects its for- 
mation history - from a two-stage binary merger, with little angular 
momentum. Box orbits do not have any net angular momentum and 
hence their 3 fundamental frequencies are (in general) uncorrelated 
with Qx < ^ ^z- This means that box orbits are generally 
found scattered in the frequency map al Q,y/Qz < 1 and to the 
left of the diagonal with ^x/^iy < 1- The long-axis (x) tubes have 
Qy ~ Qz and hence primarily cluster along a horizontal line at 
The inner long-axis tubes lie at smaller ^Ix/Clz val- 
ues near the label "J" (between the two vertical lines correspond- 
ing to the resonances (2,0,-1) and (3,0,-2)) and the outer long-axis 
tubes lie at larger ^ly/Q,z values near the label "0". The fact that 
this distribution function has only a tiny fraction of short axis tubes 
(2%) as determined by more rigorous orbital classification in VIO, 
is represented by the sparse distribution of red points along the di- 
agonal Qx/^y ~ 1. Thin dashed lines mark several possible res- 
onance lines, but we see that only the two "banana orbit" families 
labeled (2,-1, 0) and (2, 0,-1) are prominent. The qualitative dom- 
inance of the box-orbits relative to the long-axis tubes and short- 
axis tubes can be visually assessed directly from the frequency 
map, without specifically going throu gh the process of orbit clas- 
sification as done by VIO and others jCarpintero & Aguiladl 19981 : 



Long axis tubes and short axis tubes generally are not considered "reso- 
nant" orbits, however t hey may be viewed as pert urbations of their respec- 
tive "thin shell" orbits Ide Zeeuw & Huntejl990l) . 
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Figure 4. Density contours of halo mass distribution in triaxial halos after disk (shown by red contours) is grown with spin axis aligned in different ways. 
From Left to right: models SAl, LAI, lAl and TAl. Top row: density contours of halos with disk seen edge-on. Bottom row: density contours of halos in the 
plane of the disk. Halo triaxiality varies with distance from the center of the potential. 




Figure 5. Cartesian frequency maps of ~ 10'* halo particles with Vg < 200 kpc in the original triaxial halo (left); frequency map of the same pailicles after 
the growth of a disk perpendicular to the short-axis (model SAl) (middle). Right: Cartesian frequency map for particles with Rs < 10 kpc from the "sun" and 
''apo < 50 kpc. In all panels, dashed lines mark important resonances, also labeled by resonance integers (I, m, n); approximate locations of inner (3) and 
outer (0) long-axis tubes along a horizontal line Qy/Qz ~ 1 are marked. Short-axis tubes are along the diagonal line with Q^/^y ~ 1- Particles are color 
coded by energy (see text). 
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iHoffman et aTlbOOgl , l20ld) . We will see that this becomes partic- 
ularly important for identifying orbits in halos whose shapes vary 
with radius. 

Figure |5] (middle) shows the frequency map of the same set 
of halo orbits after the growth of a stellar disk perpendicular to the 
short axis of the halo (model SAl). VIO found that the fraction 
of box orbits has now decreased from 86% in the original triax- 
ial halo to 48% after the growth of the disk, but the long-axis tube 
fraction remained almost the same (12%). They found a significant 
increase in the fraction of short-axis tubes (2% to 31%) and chaotic 
orbits (1% to 9%). The increased fraction of short-axis tubes rela- 
tive to the long-axis tubes and box orbits is obvious in the frequency 
map, where the short-axis tubes appear as the enhanced clustering 
along the diagonal line. The map shows that the disk also causes or- 
bits to become separated into distinct bands in energy (color), with 
the most tightly bound (blue) particles migrating to the bottom left 
hand corner. 

The migration of particles in the frequency map occurs be- 
cause the introduction of the disk potential perpendicular to the 
z-axis, increases the vertical frequency (Qz) of orbits more sig- 
nificantly than either of the other two frequencies. Since the disk is 
centrally concentrated, the increase in is particularly significant 
for orbits which are deeper in the potential (colored blue). There- 
fore both Qy/Q,z and Q,x/^z decrease, and the blue points move 
downward and to the left. The growth of a disk also increases the 
fraction of halo orbits trapped in resonances, but some resonances 
are destroyed. For instance the vertical (2,0,-1) "banana resonance" 
has significantly increased in prominence (due to trapping of orbits 
in the x — z plane, but the (2,-1,0) banana resonance in the tri- 
axial model (which lies in the x — y plane) is destroyed because 
this is the disk plane, and the presence of the disk decreases the 
degree of triaxiality. New resonances are also populated e.g. (3,0,- 
2) the "fish" resonance, and (2,1,-2) resonancqj. Since a frequency 
map represents the ratios of the frequencies and not the frequencies 
themselves, it is insensitive to the absolute value of the energy of 
individual particles and it is therefore possible to identify the global 
orbital families and resonances (i.e. those that are important over a 
large range of energies). 

Figure |5] (right) shows a frequency map of ~ 10^ inner halo 
particles. A comparison of the right and middle panels of Figure |5] 
shows that the main features of the frequency map of the entire 
halo (middle) are also seen in the map of inner halo orbits (right). 
Thus, although the halo is more oblate within 50 kpc that at larger 
radii, the orbits in the inner halo share the major orbit families and 
resonances of the entire halo. 

The reader may be concerned that since the halo is moder- 
ately triaxial and not axisymmetric, by selecting orbits within a 
"torus" of radius Rs rather than in a region localized around the 
"sun", the frequency maps of stars in the entire torus would not re- 

^ 3-dimensi onal images of al l the m ajor resonances in triaxial potentials 
are shown in lMemtt & Valluril h999t) . 



fleet variations in the phase space structure at different locations in 
the equatorial plane. We therefore also produced frequency maps 
of subsets of the full sample of 10^ orbits, which contained only 
those orbits with initial positions within an individual quadrant of 
the galaxy model, with quadrants symmetrically about either the 
major or minor axes. These frequency maps are not shown since 
they are virtually indistinguishable from those of the full samples 
of 10** orbits (shown in the middle and right panels) demonstrating, 
that within the "solar" circle, the halo is sufficiently oblate that the 
orbital populations do not depend significantly on azimuthal loca- 
tion. 



3.2.2 LAI: Triaxial halo with long axis disk 

It has been conjectured that instead of the MW halo being oblate 
with short-axis co-aligned with the spin axis of the disk (as in 
model SAl), the disk may in fact be prolate and perpendicular to 
the long axis of the halo. Helmil ^20041 concluded, from modeling 
the tidal disruption of Sagittarius dwarf satellite, that the kinemat- 
ics of stars forming the leading arm of the Sgr stream suggest that 
the dark matter halo may be prolate with an average density axis 
ratio close to 5/3 with long axis perpendicular to the disk. This 
orientation has also been suggested by studies of distribution of 
dark matter subhalos that are satellites of MW-size d dark matter 
halos analyzed in cosmological TV-body simulations jZentner et al.l 
I2OO5I) . 

We investigated this possibility in simulation LAI, where a 
disk was grown perpendicular to the long axis of the halo. A plot of 
the projected density contours of halo particles in simulation LAI, 
after the halo relaxed into a new equilibrium with the disk potential 
are shown in Figure|6](left). The density contours show that the halo 
triaxiality varies from about 0.8 to 0.4 at small radii (a; < 20kpc), 
with long axis in the disk plane at these radii but remains triaxial 
with long axis perpendicular to the disk at larger radii. 

Figure |6] (left) shows a frequency map of 10" halo orbits ran- 
domly selected to lie with rg < 200 kpc (representing the entire 
halo DF). Recall that the frequency map of original triaxial halo 
for this model is shown in Figure |5] (left). In this map a significant 
fraction of orbits deep inside the potential (blue) populate the hori- 
zontal (1:1) resonance line corresponding to the global long-axis of 
the halo. Long axis tube orbits circulate about the a::-axis in a fixed 
direction. The dramatic increase in the length and strength of this 
family of tube orbits as well as their location on the frequency map 
(relative to the fraction long-axis tube orbits in the original triaxial 
halo in Fig.[5]left) is a direct consequence of the growth of the disk 
potential. In this simulation the disk is symmetric about the long 
X axis. Therefore orbits (especially those deeper in the potential) 
experience a somewhat larger increase Q,x, than in the other two 
frequencies. The greater increase in Qx causes the blue points on 
the frequency map to migrate towards the right of the map. Not all 
the migration to the right is associated with long-axis tubes, we see 
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Figure 6. Left: Frequency map of 10'* halo orbits in LAI with Vg < 200kpc (the frequency map of these orbits in the oiiginal triaxial halo is shown in 
Fig-ID- Axes X, y, z chosen to be the global long, intermediate, and short-axes of the halo. Right: frequency map of ^ 10* halo orbits with Rs < 10 kpc and 
apocenters less than 50 kpc. 



a significant increase in the fraction of box orbits (below the hori- 
zontal line) as well as several distinct resonances in this population. 

We also note that many tightly bound (blue) points clustered 
along the vertical line corresponding to the intermediate (y) axis 
tubes. The "intermediate-axis" tube family is generally expected to 
be unstable. However the strong clustering along the vertical line 
shows that this family is both stable and well populated in this po- 
tential. The reason for this is seen in Figure |4] (bottom row in 2nd 
column from left), which shows a slight elongation of the dark mat- 
ter (black contours) density along the y-axis for a; < 10 kpc. It ap- 
pears that in this model the j/-axis tubes appear for the most bound 
population because this is a local long axis\ 

In model LAI there is a large increase in the fraction of tightly 
bound (blue) orbits associated with the s-axis tube family, which 
happens to coincide with the symmetry axis of the disk. This is in 
sharp contrast with Figure [5] where the disk, oriented along the z- 
axis caused an increase in the fraction of tightly bound orbits asso- 
ciated with the (z) short-axis-tube family. There are also short-axis 
tube orbits in LAI but most are weakly bound (red). Thus although 
the initial triaxial halos were identical in shape and DF, the differ- 
ences in the orientation of their disk relative to the halo resulted in 
very different orbit populations (i.e. DFs) especially for the more 
tightly bound particles. 

In Figure |6] (right) we show orbits in the inner halo (Rs < 
10 kpc, Tapo < 50 kpc). The frequency map of the inner halo 
(right) is very similar to that of the outer halo (left), with the main 
difference being a decrease in the fraction of box orbits between 
the resonance lines (2,0,-1) and (3,0,-2). The colors corresponding 
to the energies of particles change purely because the range of en- 
ergies in the right hand plot is smaller. The box orbit resonances 
(below the long-axis tubes and to the right of the y-axis tubes) are 



seen more clearly, because this region is more densely populated 
with orbits due to the orbit selection criterion. 

In both SAl and LAI, the inner halos are significantly flat- 
tened in the plane of the disk (although they remain triaxial) at 
r < 50 kpc. It is therefore remarkable that the frequency maps 
of orbits confined to the inner halos of these two potentials (right 
panels of Fig. |5] and Fig. [SJ, reflect the differences in large scale 
orientation of the halos relative to the orientations of their disks. 
Although there are clearly differences between the maps of the in- 
ner halos and the full halos, these differences are significantly less 
than the differences that arise from the different orientations of the 
disk. Recall that the initial halos were identical prior to the growth 
of the disks. This implies that when accurate phase space coordi- 
nates for stars are obtained from future survey such as Gaia, they 
can be used to infer the global shape of the halo relative to the disk, 
even if accurate coordinates are only obtained within 10 kpc of the 
sun. 



3.2.3 TAl: Triaxial halo with tilted disk 

We now consider a model in which the disk was grown inclined 
to the X — y plane of the triaxial halo (model TAl). Such an ori- 
entation is motivated by studies of dark matter halos from cos- 
mological A^-body simulations that show that the relative orien- 
tation of the angular momentum axis of triaxial halos is on aver- 
age n o more than 25° — 30° from the short axis of the triaxial 
halo teailin & Steinmetjliool) . Furthermore, numerous simula- 
tions show that the disk can be misalig ned with the symmetry axis 
of the halo jvan den Bosch et al. I2OO2). 

Figure |7] (left) shows the frequency map of 10* halo parti- 
cles with Tg < 200 kpc. This frequency map is very similar to the 
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Figure 7. Left: Frequency map of 10* halo orbits with r'g < 200 kpc after the disk grows tilted at an angle of 30° to the x — y plane of the triaxial 
potential (model TAl) (the frequency map of these orbits in the original triaxial halo is shown in Fig.[5]left). Right: about 10* orbits selected from the "solar" 
neighborhood and confined to the inner halo of the model. 



model with the short-axis disk (Fig. |5] (middle)), but also shows a 
clustering of (blue) points along a vertical line al^lx/^z ~ 1 (cor- 
responding to j/-axis tubes). This family of orbits arises because the 
disk axis in this model is inclined to the x — y plane such that the 
a;-axis is in the plane of the disk but the y-axis is at an angle of 30° 
to the disk. The disk potential therefore induces resonant trapping 
of tightly bound halo orbits causing a larger circulation in the plane 
of the disk, and consequently more angular momentum about the 
y-axis. Although this family is rotating about an axis inclined to 
the y axis, this trapping appears on the Cartesian frequency map as 
an enhanced clustering of tightly bound orbits about this axis. This 
family of intermediate axis tubes does not appear in the most tightly 
bound orbits (blue) of any other model. Unlike in model LAI, the 
density contour plots do not show a noticeable elongation along the 
j/-axis. 

The DF of the inner halo of TAl (Fig.Qright) displays many 
of the features of the full DF within 200 kpc (middle), with some 
differences: the inner long-axis tube family is less prominent, the 
banana family (2,0, -1) is very sparsely populated in the inner re- 
gion and a new resonance family is seen at Q.x/^z ~ 0.6. Apart 
from these differences, the important major orbit families (the long- 
axis tubes, short-axis tubes and intermediate-axis tubes) are well 
represented, showing that the DF of the inner halo, while different 
from that of the entire halo, shares the most important characteris- 
tics that distinguish it from other models. 

Figure [3] shows that the shape of the inner halo of TAl (grey 
curves) is very similar within 50 kpc, to the shape of the inner 
halo of SAl (black curves). However the frequency map of the in- 
ner halo of TAl (Fig. [7] right) is quite different from that of SAl 
(Fig.|5]right). This implies that their DFs are very different - again 
a consequence of the different original orientation of the disks rela- 
tive to their globally triaxial halos. Thus, our analysis of the orbital 



phase space structure using the frequency maps allows us to gain 
insights into the differences in the orbital structures of two halos 
with very similar inner halo shapes. 



3.2.4 lAl: Triaxial halo with intermediate axis disk 

Recently, Law and collaborators fcaw et alj l2009l : 
iLaw & Maiewskil l2010h remeasured the shape of the MW 
halo by fitting both the velocities and positions of stars in the 
Sagittarius tidal stream with a triaxial halo model. They find that a 
slightly triaxial halo with the sun located roughly along the minor 
axis gives the best fit to the available kinematic and positional data 
of stream stars. The best-fit configuration requires the MW disk 
to be perpendicular to the intermediate axis of the triaxial halo 
- a disk configu ration that is believed to be inherently unstable 
dHeiligman & Sch warzschilj l 19791) . 

Our final model simulates this configuration with a rigid disk 
potential grown perpendicular to the intermediate axis of the triax- 
ial halo (models investigating the stability of halos with live disks 
will be presented in Debattista et al. ( 20]_1J). The frequency map of 
10* randomly selected halo orbits with Vg < 200 kpc is shown in 
Figure [H (left). The map shows a prominent clustering of weakly 
bound orbits along the vertical line corresponding to intermediate 
axis tubes with Q,x/^z = 1 at values of Q,y/Q,z > 1. Since the 
disk is symmetric about the j/-axis, at small radii the halo does be- 
come more oblate, with y as its symmetry axis. However we see 
that even the weakly bound particles (red), which are at large radii 
and expected to be less affected by the disk, are also associated with 
the intermediate (j/)-axis tube family. Figure [3] shows that although 
the shape of halo SAl (black curves) is more oblate than halo lAl 
(red curves), the radial profiles (of b/a and c/a) are very similar for 
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Figure 8. Left: Frequency map of 10'* halo orbits with Tg < 200 kpc after the disk grows perpendicular to the intermediate {y) axis in the triaxial potential 
(model lAl). The frequency map of these orbits in the original triaxial halo is shown in Figure [5](left). Right: about 10"' orbits from the inner halo of the same 
model. The j/-axis tube family is populated by weakly bound (red) points and is not important in the solar neighborhood map on the right. 



these two halos (except in the absolute degree of flattening). Con- 
sequently, when we plot only those orbits which are confined to the 
inner halo (Fig.[8]right), we see that the frequency maps of lAl and 
SAl (Fig. [5] right) are so similar that they are hard to distinguish 
from each other. The intermediate axis tube family that was seen in 
Figure [8] (left) has completely disappeared, showing that all the or- 
bits that made up this family were part of the outer halo. However, 
this similarity in the distribution functions is not entirely surprising 
since Figure|4]show that SAl and lAl have similar density contour 
distributions (this similarity was also found by DOS). 



3.2.5 Discussion of disk-halo orientation effects 

Figures[3]and|4]showed that the growth of a disk galaxy in a triaxial 
dark matter halo of arbitrary orientation modifies the shape of the 
inner part of halo, but leaves the outer part largely unaffected. How- 
ever most methods for determining the shape of the halo, assume 
that dark matter is stratified on concentric similar ellipsoids of con- 
stant shape (i.e. not varying with radius). With this assumption the 
shape of the halo can be measured with an accuracy of a few per - 
cent, out to 50 - 70 kpc jjohnston et ai]|l996l : lGnedin et al.ll2005l) . 
Since the shape of the halo probably varies significantly with radius 
due to disk formation, the assumption of a constant halo shape is 
not valid. 

Although the radial variation of the shape of the halo will be 
impossible to measure, the fact that in all cases the inner halo is 
nearly oblate and flattened like the disk will enable us to use halo 
orbits to constrain both the shape and the DF of the inner halo. 

Furthermore, the analysis of the DFs of four different halo 
models showed that a frequency map provides detailed information 
on the various orbit families that constitute the halo DF and rela- 



tive abundance of the different orbit families at various energies. In 
addition, while the disk potential traps tube orbits, which share the 
disk symmetry axis. Even restricted maps of halo orbits confined to 
the inner region (rapo ~ 50 kpc) show all the orbit families present 
in the global halo. 

It is worth noting that frequency maps are a superior method 
for identifying important orbit families in a potential than the stan- 
dard methods that rely on the properties of specific orbit fami- 
lies |CarDintero & AguilMlll99a Iwiuri et al.ll201ol ; iDeibel et all 
I2OII ). This is because the standard methods of orbit classifica- 
tion, specify the set of orbital types that will be used to classify 
orbits. They assume a priori that the shape of the halo is con- 
stant with radius. In triaxial potentials with constant shape, in- 
termediate axis tubes are unstable and not expected in the DF 
l lHeiligman & Schwarzschiljl979l) . yet we see from the frequency 
maps that three of four models contained members of this family, 
because the triaxiality of the halo varies with radius. To build a 
DF for the stellar halo, it is important to have a full representation 
of all the orbit types that are important in the halo, since the orbit 
populations, in turn, reflect the large scale shape, orientation and 
formation history of the Galactic halo. 



We caution that in all the simulations above, the disk parti- 
cles were held fixed and did not dynamically respond to the change 
in the halo. Theoretical ar guments indicate that at least som e of 
these orientations, e.g. lAl jHeiligman & Schwarzschil jl979h . are 
likely to be dynamically unstable. A more detailed analysis of self- 
consistent dynamical models of such systems are needed to as- 
certain whether such disk orientations would be found in nature 
jPebattistaetalJlOllI) . 
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Figure 9. Left: Frequency map of 10'* halo orbits with r'g < 200 kpc after a short axis disk is grown in the triaxial halo (model SA2) before the formation of 
a bar. The frequency map of these orbits in the original triaxial halo is shown in Fig.|5](left). The disk of live particles rapidly forms a bar, which subsequently 
dissolves. Middle: frequency map of the same halo orbits after the bar has dissolved. Right: frequency map of 10* orbits in the inner halo after the bar has 
dissolved. 



3.3 Non-axisymmetric halos with live disks 

Most stellar disks, including that of the MW, are not axisymmet- 
ric or time-independent since they contain features such as spirals 
and bars. These non-axisymmetric features drive angular momen- 
tum exchanges. The angular momentum exchange that occurs be- 
tween a rotating bar and the dark matter (and presumably also the 
stellar halo) is expected to result in a change in the angular mo- 
mentum of halo particles that interact with the bar ( IWeinberjl98l : 
IPeb attista & S ellwoojl20 00V In this section we study the effects 
of angular momentum and energy transfer resulting from coherent 
time-dependent perturbations from a bar. We also studied two mod- 
els in which a stellar and gaseous disk forms via realistic processes 
in a spherical halo and in a prolate halo. In these simulations, hot 
gas in the initial halo which has angular momentum, cools to form 
a rotating disk followed by star formation and feedback (see § 2 for 
details). 



3.3.1 Triaxial halo with barred disk 

To study the effects of a bar, we adiabatically grew a stellar disk that 
had only 30% of the mass of the disk in the previous simulations 
inside a triaxial halo (with symmetry axis parallel to the short axis 
of the halo, model SA2). Because of its low mass, the bar that forms 
in this syste m is eventually destroyed by its interaction with the 
triaxial halo j Serentzen & ShlosmaiJlioog) . Figure [9] (left), shows 
the frequency map of lO* halo particles with Tg < 200 kpc after 
the disk had reached full mass. The frequency map shows the two 
long-axis tube families (marked with script "J" and "0"), a strong 
short-axis tube family and several box-orbit resonances (e.g. (2,0,- 



I), (2,1,-2), (3,-2,0), as well as a few other unlabeled resonances 
seen mostly in red points). 

After the disk reached its full mass the particles were made 
"live", i.e. the disk was allowed to evolve self-consistently along 
with the halo. The rotating disk rapidly formed a bar and spi- 
ral patterns. The bar survived for ~10 Gyr and finally dissolved. 
Note that we always integrate orbits in a potential without a bar. 
This is because we "freeze" the potential before integrating orbits, 
and the presence of a bar would result in an unrealistic "freez- 
ing" of a time-dependent rotating bar. Most model s that seek to 
obtain the DF of the Milky Way disk and halo (e.g. lBinnevlbOlOl : 
iBinnev & McMillai]l201lh . neglect the bar since it is difficult to 
model its time-dependent potential. 

Figure |9] (middle), shows a frequency map of the same set 
of 10^ halo particles with rg < 200 kpc, while the right panel 
shows 10* halo particles in the inner halo, after the growth and dis- 
solution of the bar. This frequency map shows that tightly bound 
(blue) particles that were associated with the (2,0,-1) resonance 
and the short-axis tube family are scattered to other parts of the 
map. This is because these orbits are most strongly affected by the 
central bar. T he bar exchanges angular momentum with resonant 
halo particles ([l7emaine & Weinberg|l984l : lDebattista & SellwoodI 
l2000l : lAthanassoulal l2002l). Thus the blue (tightlv bound) halo par- 
ticles in the bottom right of the left panel are scattered and the res- 
onances at the bottom left of the map are no longer seen. Many of 
the blue points now appear to be associated with the outer long- 
axis tube family (at Qx/^z ^ 0.9) and the short-axis tube fam- 
ily. Resonances associated with less bound (red) particles have be- 
come broader or disappeared. Other resonance lines populated by 
the weakly bound (red) particles are only slightly affected (e.g. 2,0,- 
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Figure 10. Density contours of disk that forms in tlie spherical halo from 
cooling halo gas. The disk develops a mild oval distortion {x and y axes 
coordinates are in kpc). 

1 and 3,-2,0). These tightly bound particles are seen more clearly 
in the right panel which shows orbits from the inner halo. 

It is clear that although the middle and right plots have fewer 
strong resonances than the halo prior to the bar (left), the major or- 
bit families still appear with the same relative strengths. Although 
some resonances appear broader and those associated with parti- 
cles deeper in the potential have been scattered to other locations 
on the map, the overall structure of the frequency maps are sim- 
ilar to SAl which is also a model with a disk rotating about the 
short-axis of triaxial halo. This shows that although time depen- 
dent perturbations from a bar can affect the shape and DF of the 
halo, many resonances survive, especially those that are populated 
by orbits that do not interact strongly with the bar, and the overall 
structure of the frequency map is characterized by the global nature 
of the potential. 

3.3.2 Halos with live disks of gas and stars 

We now study frequency maps of halo orbits from models where 
a stellar disk forms from hot gas with some initial angular mo- 
mentum (defined to be about the z axis) in a spherical halo (model 
SNFWgs) and in a prolate halo (model SBgs) (see §[2jT]for details). 
Following the formation of the disk, the dark matter halo became 
fairly oblate in model SNFWgs, while the disk develops only a mild 
oval distortion (see FigllOt. In these simulations no stars form in 
the halo and therefore we once again study the orbits of dark matter 
particles. 

Figure[TT](left panels) shows the frequency maps in Cartesian 
coordinates of 10* halo stars selected with < 200 kpc in the 
initial spherical halo (top) and in the initial prolate halo (bottom). 
As we saw in several previous models (e.g. SAl, TAl), the growth 



of the disk results in an increase in the fraction of short-axis tubes 
deep in the potential (blue points along the diagonal). Both mod- 
els show prominent short-axis tube families which manifest as the 
strong clustering along the diagonal. 

VIO showed that prolate halos are dominated by long-axis 
tubes which persist even after the growth of a baryonic component, 
but these long-axis tube orbits become "rounder". Long-axis tubes 
are seen in the prolate halo as the clustering of points along the 
horizontal line with Q,y/Q.z = 1 (bottom left panel). Surprisingly 
the spherical halo also shows long-axis tubes (predominantly in the 
tightly bound orbits colored blue). This is likely to be because the 
disk forms a slight oval distortion (see Fig.llOll. Compared to the 
originally triaxial models SAl, LAI, lAI and TAI, these two mod- 
els have only a small fraction of box orbits (scattered points in the 
middle of the map) because these models are not triaxial. 

It is also instructive to analyze the same set of halo orbits in 
cylindrical coordinates (with z as the symmetry axisfl- The two 
right panels of Figure [TT] show frequency maps in cylindrical co- 
ordinates. Several resonances previously seen in Figure |2l such as 
fiz/f^H = 0.5, 1, are seen in both the spherical halo and the prolate 
halo. In addition both show new resonances n^/fi^ = ±0.5, L 
The prolate halo (bottom right) also shows a large fraction of parti- 
cles associated with the vertical resonance line at Q,^/Q,r = ±0.5, 
which are long-axis tubes. 

The disk formation process in this subsection is dynamically 
quite different from that in the previous section where rigid disks 
made up of particles were adiabatically grown in place. These two 
simulations confirm that resonant trapping of halo orbits occurs 
even during the more realistic dissipative processes by which real 
stellar disks form. 



4 USING FREQUENCY ANALYSIS TO CONSTRAIN 
POTENTIAL PARAMETERS 

In the simulations so far, we have integrated orbits in the total A'^- 
body potential for the galaxy model, which was known perfectly. 
Since the orbits were integrated in the self-consistent A^-body po- 
tential from which they were drawn, the frequency maps obtained 
represent the true DFs of these galaxy models. Since a major goal 
of current and future ga lactic surveys i s to obtain both the potentia l 
and DF of the Galaxy (lB innevll20I(j| : iBinnev & McMiIIaj|201ll) . 
our frequency based method is a promising approach for doing this. 
However, in reality the potential of the Galaxy is not known and it 
will also be measured from the spatial and kinematic distribution 
of stars. Ideally both the self-consistent DF and the potential will 
be re covered from six phase space coordinates for large numbers of 
stars teinnevll20ld:lBinnev & McMillanll20Ilh . 

^ This coordinate system is not the proper one in which to examine long- 
axis tubes, which require cylindrical coordinates where x is the axis of sym- 
metry. 
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Figure 11. Frequency maps of lO** halo particles with rg < 200 kpc in two halos in which hot halo gas cools into a star-forming disk. Left: frequency maps 
of halo orbits in Cartesian coordinates; Right: frequency maps of the same particles in cyHndiical coordinates. Top row: maps of halo orbits in the case when 
the initial halo was spheiical NFW, the bottom row shows frequency maps in the case when the initial halo was prolate axisymmetric. Both Cartesian maps 
show short-axis tubes and long-axis tubes, and only a few box orbits. Cylindrical maps show several resonances as indicated by thin dashed Hnes and labels. 



The most promising methods for measuring t he potential of 
the h a lo use kinematics of stars in tidal streams jjohnston et alj 
1 19991 : iMcMilla n & Binnev '2008*), or accurate orbits of hyper- 
velocity stars (Gnedinetal. 2005). While these methods are ex- 
pected to yield excellent estimates of the shape and density profile 
of the MW halo if it is stratified on concentric ellipsoids, all cur- 



rent numerical experiments show that the halo's shape varies with 
radius, making the measurement of the shape of the halo at all radii 
challenging due to the absence of coherent tidal streams and or hy- 
pervelocity stars over a range of radii. 

Halo stars are however, distributed over a wide range of radii 
and current studies show that even local halo stars have enough ki- 
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netic energy to travel to large radii dCarolIo et ailboiol) . We now 
test whether it will be possible to gain information about the true 
potential and DF from the frequency analysis of a large number of 
such halo stars. Jeans Theorem states that any steady state equilib- 
rium distribution function depe nds on phase space coordi nates only 
though the integrals of motion dBinnev & Tremainel2008l) . This im- 
plies that for a DF (or a random subsample thereof) that is in self- 
consistent equilibrium with its background potential, only a small 
fraction of orbits are irregular or chaotic. If a large fraction of orbits 
are chaotic (i.e. not confined to regular tori in phase space), they are 
expected to diffuse in phase space cau sing the DF to evo lve until 
the system reaches a new equilibrium jMerritt & Valluri| [l996). If 
the degree of chaoticity is low (either there are only a small number 
of strongly chaotic orbits, or there are many weakly chaotic orbits) 
the evolution of the potential co uld take some time an d fairly long- 
lived quasi-equilibria can exist jPoon & Merrittl2004l ). However, if 
the initial conditions of a large number of orbits are strongly out 
of equilibrium with the background potential, the overall degree of 
chaoticity (or "chaotic mo mentum") is large and t he self-consistent 
system will evolve rapidly jKalapotharakoslfaosI) . 

In our experiments each halo particle is treated as a test parti- 
cle which is integrated in a frozen background potential, and hence 
the orbits do not self-consistently influence the frozen potential. 
Nonetheless, the fact that we use a large (1 — 2 x 10*) ensemble 
that is a random sampling of the halo DF gives us additional col- 
lective power. A DF that is not in self-consistent equilibrium with 
the background potential w ill relax. This relax ation manifests as 
orbital diffusion or mixing jValluri et alj|2007h . This mixing will 
occur whether the potential is time independent or varying, so long 
as the DF is not in self-consistent equilibrium with the potential. 
We hypothesize that the diffusion rates log(A/) can be used to 
measure the rate of diffusion of mixing of orbits in ensemble that 
is evolved in a selected potential. If the DF is out of equilibrium 
we should obtain larger diffusion rates. We test this hypothesis by 
comparing the distributions of the diffusion rates of ensembles of 
orbits evolved in different potentials. 

Our objective is to distinguish quantitatively (with an assigned 
statistical confidence) between the correct potential and incorrect 
potentials using halo stars for which all 6 phase space coordinates 
are available. We computed diffusion rates for all 10'' orbits se- 
lected from 3 models SAl, LAI and lAl, when the orbits were in- 
tegrated in the correct potential. We also integrated the orbits from 
the DF of SAl in two "slightly incorrect" potentials: SAl rotated 
about the z axis by 10° (model SAl- 10°) and by 90° (model SAl- 
90°). Two other "slightly incorrect" potentials consisted of 10* or- 
bits from drawn from the DF of models LAI and lAl and integrated 
in the potential for SAl (models LAI-in-SAl and lAI-in-SAl re- 
spectively). 

We also considered one "strongly incorrect" model: the orbits 
of the initial triaxial halo A, were integrated in a potential consist- 
ing of halo A + a short-axis disk. Since the halo was not allowed to 
relax in response to the presence of the disk, the orbits are strongly 
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Figure 12. Cumulative distribution functions of orbital diffusion rates 
log( A/) for ensembles of 10* orbits which were evolved in the correct po- 
tential (blue curves), evolved in slightly incorrect potentials (orange curves) 
and in a strongly incorrect potential (red). 

out of equilibrium since,in effect, the disk potential is "turned on 
suddenly". 

Figure [T2lshows the cumulative distribution functions (CDFs) 
of the diffusion parameters for ensembles 10* orbits from the 3 
DFs from SAl, LAI and lAl evolved in their own potentials as 
blue curves. The CDFs of diffusion rates for 10* orbits in each of 
the four "slightly incorrect" potentials, SAI-IO°, SAI-90°, LAI- 
in-SAl and lAl-in-SAI (orange curves), and the "significantly in- 
correct" potential (red curve). The clear separation of the curves 
shows that the CDFs of log(A/) for ensembles integrated in the 
correct potential (blue curves) are always to the left of the orange 
curves implying that in the correct potentials there are significantly 
more orbits with low diffusion rates (log(A/) < —2) than when 
these ensembles are evolved in incorrect potentials (orange and red 
curves). 

We carried out pair-wise comparisons of the 8 CDFs in Fig- 
ure [12] using the non-parametric Kolomogorov-Smimov (KS) test 
to evaluate the probability p that the distributions are statistically 
different from each other. We compute the KS-statistic which mea- 
sures the "distance" d between the two distributions. Small dis- 
tances d between pairs of distributions and low (or zero) values 
of p, indicate a low probability that the two distributions are (sta- 
tistically) identical. Results of representative subsets of the tests 
are reported in Table |2] The first three tests in the table compare 
the CDF of orbits from the SAl DF evolved in the correct poten- 
tial with the same orbits evolved in three incorrect potentials. The 
small values of d and p indicate a very low probability that the blue 
curve for SAl and the orange (or red curves) are identical. The next 
two tests compare the DFs of LAI and lAl integrated in the correct 
potential (blue curves) with the CDFs when the same ensembles 
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Models compared in KS-Test 


Distance 


d 


Probability p 


SAl vs. SAl-10 


0.11 




0. 




SAl vs. LAl-in-SAl 


0.10 




4.2 X 10- 


-44 


SAl vs. Halo A + DISK 


0.15 




0. 




LAI vs. LAl-in-SAl 


9.9 X 10- 


-2 


8.6 X 10- 


-43 


lAl vs. lAl-in-SAl 


8.6 X 10- 


-2 


2.7 X 10- 


-26 


SAl-90 vs. SAl-10 


1.2 X 10" 


-2 


0.45 




SAl-10 vs. LAl-in-SAl 


1.7 X 10- 


-2 


0.12 




LAl-in-SAl vs. lAl-in-SAl 


1.9 X 10- 


-2 


0.12 




SAl-90 vs. Halo A -I- DISK 


8.8 X 10- 


-2 


1.7 X 10- 


-34 


LAl-in-SAl vs. Halo A H- DISK 


8.1 X 10- 


-2 


1.9 X 10- 


-29 



Table 2. Results of Kolomogorov-Smimov tests comparing cumulative dis- 
tribution functions of log(A/) from Figure [T2l pairvyise. 



are evolved in the potential for SAl (orange curves). These two 
tests also yielded small values of p implying a low probability that 
the CDFs of LAI and lAl were identical to those of LAl-in-SAl 
and lAl-in-SAl respectively. Thus the first 5 KS tests shows that 
the probability of CDFs of log (A/) arising from the correct and 
incorrect potentials being identical is very low (p ~ 0). 

The next 3 tests compared various "slightly-incorrect distribu- 
tions" (orange curves in Figure [T2J with each other. The results of 
the KS-tests in Table |2] show that it is difficult to statistically dis- 
criminate between the various slightly incorrect potentials because 
the probabilities of their being identical are quite large p > 0.1 

The last two tests compared two "slightly-incorrect" models 
(orange curves) with the "strongly-incorrect" models (red curve) 
(labeled "Halo A -I- DISK"): the small values of p and d indicate 
that the red curve can be statistically distinguished from the orange 
curves with high confidence. 

The tests above shows that although each particle is treated as 
a "test-particle" in the fixed background potential we are able to 
harness the collective behavior of an ensemble drawn from a self- 
consistent distribution function assumed to be in steady state to 
statistically identify cases where the orbit ensemble is not in equi- 
librium with the background potential. Although the diffusion rate 
of an ensemble (in a fixed background potential) can be quantified 
even without self-consistent calculations, it is more correct to think 
of this diffusion in terms of the collective relaxation that occurs 
when an N-body system is out of equilibrium, rather than chaos, 
since these orbits are drawn from a distribution function. 

This is a fundamental result that relies on the Jeans Theo- 
rem: when the initial positions and velocities (phase space coordi- 
nates) of orbits are not drawn from a self-consistent DF in a steady- 
state equilibrium potential, most orbits are not launched on regular 
tori and hence they will diffuse in phase space (trave l through the 
"Arnold web", c.f. iLichtenberg & Liebermanlll992h . The results 
above show that although there is a spread of diffusion rates in any 
given ensemble. Since the CDF of the distribution for the entire en- 
semble is shifted to smaller values of the diffusion rate when the 



potential is close to equilibrium and becomes statistically larger 
when the potential is incorrect. 

From Figure [12] we see that the CDF of diffusion rates of a 
"slightly incorrect" model can be clearly distinguished from the 
correct model. This suggests a novel way to utilize the six phase 
space coordinates of MW halo stars to distinguish between various 
possible models for the halo potential. However, the overlap of or- 
ange curves for various "slightly incorrect" models indicates that it 
may not be easy to distinguish between the various "slightly incor- 
rect" models in a way that allows us to progress iteratively towards 
a model closest to the true potential. 

Nonetheless this method has the value that it will work best in 
the region of the MW halo that is more well mixed (i.e. the inner 
halo), whereas methods such as modeling the tidal tails of dwarf 
satellites relies on their being less well mixed and will primarily be 
applicable in the outer halo. Note that it is not necessary to assume 
that all the orbits in the correct equilibrium model are regular. In 
fact, the CDF of diffusion rates in Figure[72lshows that even in the 
correct/equilibrium models there are a few orbits which have fairly 
high diffusion rates. 

It is important to note that this method will not work if the net 
potential is spherical (which the MW is not because ~10% of the 
total mass is in the highly flattened Galactic disk), or if the potential 
is assumed to be of Stackel form. This is because all orbits in spher- 
ical and Stackel potentials are integrable by definition (i.e. all or- 
bits are perfectly regular), hence the diffusion of orbits would only 
measure numerical errors of the method. However, since Stackel 
potentials are idealized potentials of primarily theoretical interest, 
there is no reason to expect that the Galaxy should have such a 
form. In a future paper we will refine these ideas to quantitatively 
assess the possibility of measuring the shape of the Milky Way's 
inner halo with ensembles of stellar orbits. 



5 DISCUSSION AND CONCLUSIONS 

Understanding the structure, dynamics and formation history 
of the MW is a major thrust of current astronomical research. 
The fundamental science driver of the burg eoning Gala c tic all - 
sky survey industry e.g. SDSS - SEGU E jYannv et alj l2009h, 
APOGEE jAUende Prieto et al.' 2008"), RAVE (ISteirmietzl 
r2003), LSST ILSST Science Collaborations et alj |2009|), 



PanSTARRS jKaiser et alj 2002), LAMOST (Hu& Jiang 200l^ 
Skymapper jKeller et alfBoO?.). HERME S (iBarden et al., ,2001.. 
and eventually Gaia jPerrvmanet"ai]|200lh ) is the premise that the 
Galaxy and the Local Group were formed in a manner that typifies 
the formation of galaxies in the ACDM paradigm. Therefore 
understanding the structure, dynamics, chemistry and thereby 
formation history of our own Galaxy will result in stronger 
constraints on the complex physics of galaxy formation. 

The use of orbital frequencies of halo stars has been recog- 
nized as an important way to identify t he relics of satellites which 
were tidally disrupted in the MW halo jMcMillan & Binnevll20o3 : 
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iGomez & Helirillioioh . In a recent paper iGomezet al] i2010l) also 
showed that when the time dependent growth of the Galaxy and 
observational errors associated with Gaia are taken into account, 
it will be possible to uniquely identify about 30% of the accre- 
tion events that occurred in the last 10 Gyrs, while the remainder 
of the accretion events will likely be difficult to disenta ngle from 
a more smoothly mixed component. iHelmi et"aL ll) have ar- 
gued, from a comparison of the degree of spatial variation in the 
distr ibution of stars i n the stellar halo observed in the SDSS-II sur- 
vey teell et al.ll2003) and the simil ar distributions of h alo stars in 
cosmological A'^-body simulations dCooper et al1l201(]|) . that there 
may exist a smooth underlying stellar halo component that is much 
more well mixed than material solely accreted from satellites. The 
method outlined here complements the work that focuses on disen- 
tangling the relics of tidal streams, in that it can be applied to both 
well mixed and unmixed orbits and can therefore be applied to a 
much larger sample of halo stellar orbits. 

The SDSS-SEGUE survey has already obtained phase-space 
coordinates for over 17,000 stars within 4 kpc of the sun. In a paper 
in preparation, Derris et al. ( 2011) apply the techniq ues described 
in thi s paper to the SDSS-Segue Calibration sample jCarollo et al.l 
l2010i) . and correlate the orbital properties of halo stars with their 
metallicities and spatial distributi on. They use revis ed distances to 
the stars in the Calibration sample (i Beers et"al ]20T 1) that overcome 
recently reported distance measurement errors ( Schoe nrich et al.l 
I2OIOI) . Upcoming surveys will increase both the volume of space 
observed and the accuracy with which distances, proper-motions, 
radial velocities and metallicities are measured for MW halo stars, 
significantly impacting our understanding of the structure and dy- 
namics of the Galaxy. If the stellar halo does consist of distinct 
components — a well-mixed inner halo which formed largely in 
situ (i.e. in potentials deeper than typical dwarf satellite potentials) 
and an outer halo that was formed primarily from the accretion of 
tidally stripped dwarf satellites — one might also expect to observe 
distinct differences in the orbital populations of these two compo- 
nents. 

One of the primary goals of current and future surveys of the 
MW is to determine the shape and radial density profile of the dark 
matter halo and to construct the self-consistent phase space dis- 
tribution function of the major stellar components - the thin and 
thick disks, the bulge and the stellar halo (Binney & McMillan 
l201lh . A popular and flexible method for constr ucting the DF 
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of a poten tial is th e orbit superposition method |Schwarzschilc 

Ivan der Mare l et al.lli998l : ICretton et al] 1 19991 : 1 Valluri et al 

2004lThomasetalJ 20041) . which relies on orbit libraries that uni- 



formly sample orbital initial condition space. In the DF of realis- 
tic triaxial potentials like the halos of galaxies, a large fraction of 
boxlike orbits are resonant (i.e. have commensurable frequencies ) 
jMiralda-Escude & Schwarzschilj[T98^ : iMerritt & Vallurilll999l) . 
Determining the fraction of orbits associated with resonances is 
therefore important for constructing the distribution function. How- 
ever, since these orbits are "resonantly trapped", they densely pop- 



ulate very small regions of phase space, and are likely to be under- 
represented in the orbit libraries used to construct DFs (which gen- 
erally sample some initial condition space uniformly). Resonant or- 
bits are also not adequately represented in the DFs constructed via 
orbital torus construction methods, since they are, by definition, not 
present in perfectly regular Stackel potentials, whose orbital tori 
are adiaba tically deformed to construct the DF of realistic poten- 
tials (e.g. Binnevll2010l :lBinnev & McMillaij |2oTTh . The methods 
described in this paper address this important issue in a uniquely 
powerful way. Rather than attempting to guess the initial distribu- 
tion of orbits required to populate the orbit libraries, the frequency 
analysis of the orbits of halo stars in a set of trial potentials can 
be used to construct the frequency map, which we have shown to 
yield a robust picture of the global DF, even when orbits are se- 
lected within a limited volume in the galaxy. The measurement of 
the diffusion rates of the orbits then gives an estimate of the best 
trial potential, and can be used as an additional constraint in the 
Schwarzschild or torus construction method. The applicability of 
these methods to the construction of the full phase space distribu- 
tion function of the stellar halo, will be discussed in greater detail 
in future work. 

In the majority of simulations presented in this paper we 
do not consider the possibility that the inner and outer halos 
could have angular momenta of different directions or signs. 
Warps seen in disk galaxies, are now believed to be evi- 
dence that the angular momentum vectors of outer dark mat- 
ter halo is misaligned f rom the disk and the associated in- 
ner halo (e.g. .Ostriker & Binnev 19891: Jiang & BinnevI 199^: 



Debattista & Sellwoodll999l : IShen & Sellwoodl2006l : lRoskar et all 
2OIOI) . These authors argue that cosmic infall causes the angular 
momentum direction of the outer dark matter halo to be different 
from that of the inner dark matter halo. The build up of dark matter 
halos via cosmic infall occurs due to the accretion of satellites, and 
such hierarchical infal l is also believed to build up the stellar halo. 
In a companion paper ( I Valluri et al.l201 lb we study the angular mo- 
mentum distributions of both star particles and dark matter particles 
in high resolution cosmol ogical simulations of disk galaxies from 
the MUGS collaboration dStinson et al.ll2010i) . Preliminary results 
indicate that the angular momentum distributions of halo stars at 
different radial distances from the galactic center can be used to 
trace the angular momentum distribution of the dark matter halos 
at these radii. Thus orbital analysis of halo stars could potentially 
be used to measure the angular momentum of the dark matter halo 
and test the prediction that differential angular momentum of the 
halo can produce warps. Model TAl (Figures H (right) &|71( show 
the effect of a disk inclined to the equatorial plane of a triaxial 
halo. Figure |4] (right) shows that the density contours of the inner 
halo become aligned with the disk, while the outer halo remains un- 
changed. In this simulation the disk was held rigid, but it is likely 
that a live disk would form a warp due to the misalignment of the 
inner and outer halo angular momenta. Nonetheless this simulation 
is illustrative of the possibility of using frequency maps of halo 
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stars to detect such tilts in the principle axes and angular momen- 
tum directions of different parts of the halo since it illustrates that 
the orbital frequencies only depend on the assumed potential. 
To conclude we summarize the main results of this paper. 

(i) Frequency analysis of the orbits of a large representative 
sample of halo particles can be used to construct frequency maps 
which provide a compact representation of the distribution function 
of the stellar or dark matter halo. The map also give an estimate of 
the fractions of orbits in different major orbit families and enables 
easy identification of dominant resonances. The identification of 
global resonances is important for constructing global DFs for the 
MW since global resonances strongly constrain the DP. The ability 
of frequency maps to reveal globally important orbit families and 
resonances, even with orbits selected from a limited volume around 
the sun, suggests that this method will provide sig nificant input to 
the co nstruction of DFs for the Milky Way galaxy jMav & BinnevI 

[mi). 

(ii) The adiabatic growth of a disk in a halo results in significant 
resonant trapping of halo particles with resonant orbits appearing 
clustered along narrow resonance lines on a frequency map. In a 
Cartesian frequency map of orbits in a moderately triaxial halo we 
see several resonant box orbit families. In axisymmetric halos, the 
use of cylindrical coordinates reveals a large number of resonances 
primarily between the radial frequency Q.11 and the vertical fre- 
quency Q,z. The strong resonant trapping seen in all cases where the 
disk grows quiescently in a pre-existing halo imply that the identi- 
fication of resonances in the Milky Way's stellar halo could provide 
evidence for the adiabatic growth of the disk following the forma- 
tion of the halo. 

(iii) Resonances are found in all of the controlled simulations 
in which disks were grown inside halos, regardless of the shape 
of the halo and regardless of how the disk was oriented relative 
to the large scale orientation of the halo. It has been previously 
demonstrated that the growth of a bayonic disk inside a triaxial 
halo deforms the inner regions to make them more oblate, but the 
halos remain modestly triaxial at large radii. We find that although 
the inner regions of the halo are similar regardless of large scale 
halo orientation, the frequency maps of inner halo particles reflect 
the differences in the global orientation of the halo relative to the 
disk. 

(iv) We see that halo resonances are formed by both static disks 
and live disks (those that form spiral features and bars). However, 
the coherent time dependent perturbations from a bar can result in 
scattering of the most tightly bound particles, resulting in broader 
(less well defined) resonances. 

(v) Controlled hydrodynamic simulations, in which hot gas dis- 
tributed throughout a dark matter halo, is allowed to cool and form 
a gas disk which then forms stars, give rise to halo frequency maps 
with resonant structure, much like the adiabatic simulations. This 
shows that resonant trapping of halo stars by the disk are not purely 
a feature of idealized coUisionless simulations in which a rigid disk 
is grown in place. 



(vi) We find that the cumulative distribution of orbital diffusion 
rates are lower by a statistically significant amount, when a large 
ensemble of orbits is integrated in a potential in which it is in self- 
consistent equilibrium, and that the orbital diffusion rates are sig- 
nificantly larger when the potential is incorrect. This can potentially 
provide a novel way to constrain the potential of the MW directly 
from the 6-phase-space coordinates of a large sample of Milky Way 
halo stars. 
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